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• Abstract 

We analyze the simultaneous time-optimal control of two-spin systems. 
\Q | The two non coupled spins which differ in the value of their chemical off- 

sets are controlled by the same magnetic fields. Using an appropriate 
rotating frame, we restrict the study to the case of opposite shifts. We 
then show that the optimal solution of the inversion problem in a rotat- 
ing frame is composed of a pulse sequence of maximum intensity and is 
similar to the optimal solution for inverting only one spin by using a non- 
resonant control field in the laboratory frame. An example is implemented 
experimentally using techniques of Nuclear Magnetic Resonance. 

cr 

1 Introduction 

. Since its discovery in 1945 by Purcell, Torrey and Pound, Nuclear Magnetic 

j^. ' Resonance (NMR) has become a powerful physical tool to study molecules and 

matter in a variety of domains extending from biology and chemistry to solid 
physics and quantum mechanics pQ. NMR involves the manipulation of nuclear 
. spins via its interaction with a magnetic field, and is therefore a domain where 

techniques of quantum control can be applied (see [2] and references therein). 
Such an approach has many potential applications ranging from the improve- 
ment of the resolution and sensitivity of NMR spectroscopy experiments [3] to 
quantum computing [4 . The control technology developed over the past fifty 
years allows the use of sophisticated control fields for spectroscopy and also 
permits the implementation of complex quantum algorithms [5J. 
' In this context, some challenging control problems are raised by the experi- 

mental constraints of NMR experiments. Roughly speaking, the measured signal 
is the magnetization of a sample which is produced by a large number of spin 
systems. One usually assumes in simple models that the static magnetic field is 
the same across the sample, i.e. the field is perfectly homogeneous with respect 
to the different spins. This is not always true in practice since for technical rea- 
sons it is difficult to generate homogeneous fields. Even in the situation where 
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the magnetic field is uniform on a macroscopic scale, the interaction between 
the different atoms (or between a spin and its environment) induces a chemical 
shift on the frequency transition of a given spin. This leads classically to an 
unwanted rotation of each individual spin around a fixed axis which is not taken 
into account in the simplest model of spin 1/2 particles. The shift being differ- 
ent for each spin, the rotation is different for each spin. Note that this effect is 
useful in NMR spectroscopy since it encodes in a sense some informations about 
the structure of the molecules. The consequences are negative from a control 
point of view since this phenomenon decreases the efficiency of the control field. 
The objective is therefore to find controls able to bring the system towards a 
given target state in a sufficiently robust way with respect to inhomogeneities of 
the transition frequency. This problem has been solved numerically in different 
works [15| leading to very efficient but complicated solutions. In particular, no 
insight into the control mechanism is gained from this approach and no opti- 
mality result has been proven. Note that some related works have been done 
in the control of molecular dynamics by laser fields |16| by using monotonically 
convergent algorithms |17j . 

In this paper, we propose to revisit this problem by using techniques of 
geometric optimal control theory [7J . Geometric optimal control is a vast 
domain based on the application of the Pontryagin Maximum Principle (PMP) 
where the idea is to use the methods of differential geometry and Hamiltonian 
dynamics to solve the optimal control problems [51 [7] . This geometric framework 
leads to a global analysis of the control problem which completes and guides 
the numerical computations. Some geometric results on the optimal control of 
spin systems have been first obtained by N. Khaneja and his co-workers |12| . 
Recently, the time-optimal control of dissipative spin 1/2 particles has been 
solved theoretically |14| and implemented experimentally [13]. In this work 
we study the simultaneous control of two non-interacting spins with different 
resonance frequencies. More precisely, we consider as an example the problem 
to simultaneously invert the magnetization vectors initially aligned along the z- 
axis defined by the direction of the static magnetic field. 

Using an appropriate rotating frame, we show that we can always consider 
the symmetric case where the two transition frequencies are opposed. In this 
situation, the time-optimal solution for inverting the two spins by the same 
transverse radiofrequency (rf) control fields is a bang-bang pulse sequence in a 
frame rotating at the rf frequency. The remarkable point is that the co-rotating 
component of the applied rf field is the same as the one used to invert only one 
spin with one non-resonant control field in the laboratory frame [9]. We finally 
implement experimentally the optimal solution by using techniques of NMR. 

The paper is organized as follows. In Sec. [51 we recall the tools to control one 
spin in minimum time with a transverse magnetic field which is not in resonance 
with the frequency of the spin. In Sec. [31 we establish that this control field is 
also the optimal solution to simultaneously invert two spins. An experimental 
illustration is given in Sec. 3) A summary of the different results obtained is 
presented in Sec. [SJ 
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2 Time-optimal control of a spin 1/2 particle 
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We consider the control of a spin 1/2 particle whose dynamics is governed by 
the Bloch equation: 

\ 

-uj x M z (1) 

UJ X My J 

where M is the magnetization vector and ui the chemical shift offset. The dy- 
namics is controlled through only one magnetic field along the x- axis which 
satisfies the constraint uj x < uj max . We introduce normalized coordinates x — 
(x,y,z) = M/Mq where Mq = M$e z is the thermal equilibrium point, a nor- 
malized control field u x — 2iru} x /u) max which satisfies the constraint u x < 2ix 
and a normalized time r = uj rnax t / (2t:) . Dividing the previous system by 
w m axMo/(2ir), we get that the evolution of the normalized coordinates is given 
by the following equations: 

I / -Ay \ ( \ 

(2) 

where A is the normalized offset given by A = 2ttuj / uj max . 

The complete description of the time-optimal control problem of a spin 1/2 
particle by a non-resonant magnetic field is done in Ref. [9J . In this section, we 
give only a brief summary of the results of this paper which will be used in our 
study. The reader is referred to [H| for the different proofs of these results. Note 
that when the spin is controlled by two magnetic fields, respectively along the x 
and y directions, then the system is equivalent to a two-level quantum system in 
the rotating wave approximation [8 J . This means that a unitary transformation 
can be used to remove the drift term depending on A. In this case, the optimal 
control field is a 7r-pulse. 

The problem we consider belongs to a general class of optimal control prob- 
lems for which powerful mathematical tools have been developed [TU]. They 
correspond to systems on a two-dimensional manifold (here the Bloch sphere) 
controlled by a single field. The evolution of the system is ruled by the following 
set of differential equations: 

f = F(x) + uG{x) (3) 

where x is the two-dimensional state vector and u the control field which satisfies 
the constraint u < uq with, here, uq = 2tt. The time-optimal control problem is 
solved by the application of the Pontryagin Maximum Principle (PMP) which 
is formulated using the pseudo-Hamiltonian 

H =p-(F + uG)+p Q , 

where p is the adjoint state and po a negative constant such that p and pa are 
not simultaneously equal to 0. The Pontryagin maximum principle states that 
the optimal trajectories are solutions of the equations 

H(x,p,v) = max|„|< Uo H(x,p,u) (4) 
H(x,p, v) = 0. 
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Introducing the switching function $ given by 



m=p-G, 

one deduces using the second equation of Q that the optimal synthesis is com- 
posed of concatenation of arcs 7+, 7_ and 7s- 7+ and 7_ are regular arcs 
corresponding respectively to sign[$(t)] = ±1 or to the control fields u — ±uq. 
A switching from uq to — uq or from — uq to uq occurs at t = to when the 
function $ takes the value zero and when this zero is isolated. Singular arcs 
7 S are characterized by the fact that $ vanishes on an interval [to,ti]. I n this 
case, differentiating two times $ with respect to time and imposing that the 
derivatives are zero, one obtains that the singular arcs are located in the set 

S = {x; A s (x) = det(G, [G,F])(x) = 0}. 

We recall that the commutator [F, G] of two vector fields F and G is defined 
by: 

[F, G] = VF • G - F • VG 

where V is the gradient of a function. The singular control field u s can be 
calculated as a feedback control, i.e. as a function of the coordinates by imposing 
that the second derivative of $ with respect to time is equal to 0: 

[G, [G,F}}+u s [F,[G,F}}=0. 

The optimal solution can follow the singular lines if the control field is admissi- 
ble, i.e. if |u s (x)| < u . 

Since the two-dimensional manifold of our control problem is the Bloch 
sphere, the adapted coordinates are the spherical ones: 

r sin 9 cos <f> 

r sin# sin0 , (5) 
r cos 6 

which leads to the following system: 

/ \ 

(6) 

v — cos 4> cot 9. I 
The pseudo-Hamiltonian H has the form: 

H = Ap^ — u(sin <j>pQ + cos <j) cot 0p$). (7) 
and the switching function is given by $ = simppg + cos <f> cot 9p^. Since 

( ° 

[G, F] = -A cos <p 

\ A sin 4> cot I 

one deduces that S is the set 

S = {x\ sin^cottf = -cos 2 (f> cot6>} = {x\9 = tt/2} , 




i.e. the singular locus is the equator of the Bloch sphere. The time-optimal 
control problem is solved in [3]. It has been shown that the optimal solution 
reaching the south pole from the north pole is the succession of different bang 
pulses, i.e. of pulses of maximum intensity 2n. The number of bangs is at most 
equal to 2 if A < 27r and can be larger if A > 2tt. The singular extremals play 
no role for this spin inversion. An example of an optimal pulse sequence and the 
corresponding trajectory is displayed in Figs. Q]and[5J The optimal trajectory 
is not smooth at switching points where the value of the control field changes. 
The switching times can analytically be determined using the material of Ref. 
[9]. In the case of Fig. [I] the optimal solution is a type-2-trajectory described by 
Proposition 5 of [9]. These times can also be computed numerically by solving 
a shooting equation. More precisely, this means that one has to determine the 
initial adjoint state p(0) ~ (p#(O),p0(O)) such that the corresponding Hamilto- 
nian trajectory (x, p) with initial conditions (x(0),p(0)) goes to the target Xf 
at time tf. This condition can be expressed as the determination of the roots 
of the equation x(tf)[p(0)} — Xf which can be solved if one has a sufficiently 
good approximation of p(0) by a standard Newton-type algorithm. Note that 
the control field is determined along the trajectory by computing the switching 
function $. 

At this point, we can extend the previous discussion as a first step towards 
the simultaneous inversion of two spins. We analyze the dynamics in a rotating 
frame by using the rotating wave approximation (RWA) where the offsets of the 
two spins are symmetric and given by ±A. The rf field is assumed to be at the 
rotating frame frequency. As a consequence of the symmetries of the problem, 
one sees that if u(t), the co-rotating component of the applied rf field, steers 
the spin with offset A from the north pole to the south pole then the same field 
will also invert the other spin. The trajectories of the two spins in the y- and z- 
directions will be the same, while it will be opposite along the x- axis. Note that 
this solution is not the unique solution and a family of solutions satisfying the 
same requirement can be determined. Consider the set of control fields defined 

by 

u x = u(t) cos a 



= u{t) sin a 



(8) 
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— sin a 
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where a € [0, 2n]. If we consider the following rotation R(a) of angle a along 
the z- axis: 

' X 

(9) 



then the new system in the coordinates (X, Y, Z) is controlled by a single field 
u(t) along the X- direction. It is also straightforward to see that this solution is 
the optimal one for the inversion control of two symmetric spins by one control 
field. 

The question that we ask now is whether this simple solution is the optimal 
solution of the simultaneous inversion of two spins when two control fields are 
considered. 
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3 Simultaneous control of the inversion of two 
spin 1/2 particles 

3.1 The model system 

We consider two different spin 1/2 particles with the offsets u> a and u>b- Using 
the same normalization as in Sec. [5] and the RWA, one arrives to the following 
equations: 
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(10) 



where the coordinates (x a ,y a , z a ) and (xb,yb,z ) are respectively associated to 
the first and second spins a and b. The parameters A a and A& are the offsets 
of the spins a and b with respect to the frequency of the rotating frame. The rf 
field is here also at the rotating frame frequency. We assume that the two spins 
have the same equilibrium point Mq. As mentioned below, two magnetic fields 
along the x- and y- directions are taken into account in this problem. They 

satisfy the constraints 



. +Uy < 2ir. Using a rotating frame that rotates at 

frequency (A + A&)/2, it is straightforward to transform this system into a 
symmetric one where the frequencies of the two spins are opposite. This is the 
case analyzed below. 

We introduce the spherical coordinates for the two spins and we get: 



n 



V A J 



( o \ 


A 



V -a ) 



+ U X G X + UyGy 



(11) 



where 



V 



o 

- sin (p a 

COt 9 a COS (f) a 



— sin <f>b 
cot 9b cos 4>b 



G y — 



\ 



\ 

COS 4> a 

cot 9 a sin (j> a 


COS <pb 

cot 6b sin <pb } 



Since the radial coordinates (r a ,rb,p ra ,Pr b ) play a trivial role in this problem, 
we omit them in the following equations. 

We apply the PMP to this system in the time-optimal case and we obtain 
the following pseudo-Hamiltonian 



H = 



P<t>b) - {U X G X +UyGy) 



(12) 



G 



where p is the adjoint vector of coordinates (pg a ,p t j >al pg b ,p l j )b ). In the normal 
case, the optimization condition leads to the following optimal controls: 

V'G X p ■ G y 



y/(p-G x y + (p-G v )*' v y/(p-G x )* + (jt-G y )*' 

where p ■ G x and p ■ G y are not simultaneously equal to 0. The singular case 
occurs when p ■ G x = p ■ G y = 0, which defines the switching surface S. In 
the two-control problems, singular trajectories are the trajectories which lie on 
S. We assume in this paper that these controls do not play any role in our 
problem. This is expected since singular extremals are not generically optimal 
for a two-control problem [TS] . 

We get the normal Hamiltonian H n by replacing the control fields by their 
expressions: 



(15) 



H n = A( P4>a - P4>b ) + G x ) 2 + (p. G v f. (14) 

The normal extremals will be given by the Hamiltonian trajectories of H n . The 
next step of our study will consists in the analysis of this Hamiltonian flow. 

For that purpose, we introduce the following canonical transformation on 
the <t>- coordinates: 

</>b 
</>b 

which is defined through the generating function: 

F 2 = t;P4>A<P+ + </>-) + -jP4> b {4>+ - <t>-) 

with the transformation: 

0F 2 dF 2 , dF 2 , dF 2 



This leads to 



P<t>+ ~ P4>- 



(16) 



The Hamiltonian H n expressed in the new set of coordinates does not depend 
of (j) + , so p^ + is a constant of the motion. Since at the initial time in the north 
pole, p<f, a {Q) = P(p b (0) = 0, one deduces that p l f )+ = 0. One finally arrives at 

H n = 2A P4> _ + [pl a + p 2 0b + pl_ (cot 2 9 a + cot 2 e b ) 
+2 cos 0_ (pe a Pe b - p\_ cot a cot 9 b ) 
-2p c /,_ sin (j>^(pe a cotd b - p 8b cot 6» a )] 1/2 . 

Care has to be taken with the use of these coordinates on the poles of the 
sphere. On a pole, we have cot 9 — > ±oo and p^ = but the product p^, cot 9 
remains finite. In this paper, spherical coordinates are only used to describe the 
geometric properties of the extremals and to highlight their symmetries. All the 
numerical computations are done in cartesian coordinates. 

Note also the symmetric role played by 9 a and 9 b in the Hamiltonian H n . 
This symmetry will be used in the proof below. 
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3.2 The optimal control problem 

We first analyze the characteristics of the extremal trajectories which are so- 
lutions of the control problem. In particular, if the inversion is realized by an 
extremal trajectory then the following relations are satisfied: 

Vf G [0,t f ],pg a (t) = pe b {t) and 6 a (t) = 6 b {t), 

where tf is the control duration. 

To show this property, we assume that the south pole is reached by the 
extremal. In this case, the final point satisfies by definition: 

Oa(tf) = 9 b (t f ), e a (t f ) = 6 b (t f ), P4> _ = 0. 

Using the Hamiltonian H n , we obtain: 

{ 6a = = (PSa + cos $-Pe h ~ P<p- sin 0_ cot 6 b )/ 



where 



Iff = (Pe b + cos <i>-Pe a ~ P<f>- sin 0_ cot 9 a )/ 



(17) 



Q=Pe a + Pe b + ?4_ (cot z 9 a + cot" 
+2 cos 0_ {pe a pe b - P%_ cot 8 a cot 8 b ) 
-2p tj) _ sin 0_ (pe a cot 9 b - pe b cot 6 a ). 

We consider that Q{0) ^ which is always possible by a judicious choice of the 
initial adjoint vector p(0). This implies that H n > since p l j ) _ — at a pole. 
Using the fact that H n is a constant of motion, one deduces at the final point 
that Q(tf) 0. From Eqs. (fTT|) , one finally arrives at 

(P8 a (tf) -Pfl* (*/))(! - cos <£_(*/)) = 0. 

If cos <fi-(tf) = 1 then Q = 0, which is not possible from our hypothesis. We 
have therefore pe a (tf) — Pe b (tf) and using the Hamiltonian equations, we then 
obtain that 9 a (t) = 8 b (t) and pg a (t) — pe b (t) for any time t since at time t = tf 
we have 9 a (t f ) = 9 b (t f ) and pe a (t f ) = pe b (t f ). 

We also get that cp + = 0, i.e. cf> + = <fi+Q. In the new coordinates X a 
and X b such that X a — i?(0 + o/2)a; a and X b — R((f> + Q/2)x b , the sum of the new 
aziinuthal angles is zero and we obtain the following symmetry on the trajectory: 

X a (t) = X b (t) 

Y a (t) = -Y b (t) (18) 
Z a (t) = Z b (t) 

for any t G [0,£/]. In these coordinates, the two control fields are given by 



u x {t) = u x (i)cos(f )+u y (t)sm{!f) 
u Y (t) = -u x {t) sin(^) +%(*)cos(^ 



(19) 



From this symmetry, we have ux(t) = and thus u x (t) cos(0o+ /2)+u J/ (£) sin(0o+ /2) 
0. Since u x {t) 2 + u y (t) 2 < Ait 2 , this leads to 

u x (t) = u (t) cos(</>o+/ 2 ) (2Q) 
u y (t) = -u (i) sin(0 o +/2) 

where uo(t) is a bang-bang pulse of amplitude 2tt. We therefore recover the case 
of Sec. [5] of the optimal control of one spin system. 
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Figure 1: (Color Online) Panel (a): Plot of the optimal control field for 
the simultaneous inversion of two spins with offsets u}/2ir = ±483 Hz and 
Umax/Zn = 120.75 Hz. Panel (b): Plots of the corresponding optimal trajecto- 
ries for M°(t) (solid black line) and My(t) (red dashed line). The experimentally 
measured trajectories are presented by open circles. Panel (c): Plot of the opti- 
mal trajectory for M"(t), together with the experimental representations (open 
circles) . 
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4 Experimental illustration 



Here we demonstrate the inversion for the case where the symmetric offsets 
tj a = —LOb are four times larger than the maximum radiofrequency (rf) amplitude 
using the techniques of nuclear magnetic resonance. The optimal pulse (a bang- 
bang pulse) is shown in Fig. Q]and implemented on a Bruker Avance 600MHz 
spectrometer with linearized amplifiers. The experiment was performed using 
the two distinct proton spin signals of methyl acetate (dissolved in deuterated 
chloroform). The two resonances, one from the — OCH3 moiety and the other 
from the — OOCCH3 moiety, were separated by 966 Hz in the 1 H NMR spectrum. 
The irradiation frequency was positioned in the center of the two peaks, i.e. 
ujq = (uj a + uji,)/2, resulting in offsets of uj = (ui a — w&)/2 = 2ir x 483 Hz for 
the two resonances. The maximum rf amplitude was chosen to be uj max = 
lo/A = 2ir x 120.75 Hz and the duration of the optimal inversion pulse shown 
in Fig. Q] is T p =6.409 ms. At room temperature (298K), the experimentally 
measured relaxation time constants of the two spins are Tf w T\ = 4.95 s, 
Tg w T| = 140 ms, which have a negligible effect during the much shorter 
pulse duration T p . The x and y components of the Bloch vectors, M x ,b (t) and 
My ,b (t) were measured experimentally by interrupting the optimal pulse shape 
after the time t and measuring the amplitude and phase of the signal after 
Fourier transformation of the resulting free induction decay (FID). In order to 
measure the z component of the Bloch vectors, the experiments were repeated 
with the addition of a pulsed magnetic field gradient (of duration about 0.2 ms 
with sine shape), followed by a 90° hard pulse. A reasonable match between 
simulated and experimentally determined trajectories is found. For example, 
Fig. [1] shows the simulated and experimental trajectories of M£(t), My{t) and 
M"(t) as a function of time. Fig. [5] shows the projections of the simulated and 
experimental trajectories of both Bloch vectors. 

5 Summary 

In this last section, we give a brief overview of the results obtained in this 
paper. The four relevant cases for the simultaneous inversion of two spins are 
the following. 

1. two control fields along the x- and y- directions and one offset uj: the 
optimal solution is a n- pulse [8]. 

2. one control field and one offset uj: the optimal solution is a bang-bang pulse 
sequence with a number of switching depending upon the ratio ui/u> max 

m- 

3. one control field and two offsets uj and —uj: The optimal solution is the 
same as in (2). 

4. two control fields and two offsets u) and — uj: The optimal solution is also 
the same as in (2). 
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Figure 2: (Color Online) Projections of the optimal trajectories for the inversion 
of the two Bloch vectors with offsets, oj/2-k = ±483 Hz, in the (x, z)- and (y, z)- 
plane are shown in (a) and (b), respectively. The simulated trajectories of the 
two spins are plotted in black solid line and red dashed line. The experimentally 
measured trajectories of the two magnetization vectors are represented by open 
and filled circles. 
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